Coarse graining the Cyclic Lotka-Volterra 
Model: SSA and local maximum likelihood 
estimation. 



C. P. Calderon 1 , G. A. Tsekouras 2 ' 3 , A. Provata 2 , and I. G. Kevrekidis 1,4 

1 Department of Chemical Engineering, Princeton University, Princeton, New 
Jersey 08544-5263, USA 

2 Institute of Physical Chemistry, National Research Center "Demokritos" , 15310 
Athens, Greece 

3 Physics Department, University of Athens, Panepistimioupolis, 10679 Athens, 
Greece 

4 Corresponding Author: yannis@arnold.Princeton.edu 

Summary. When the output of an atomistic simulation (such as the Gillespie 
stochastic simulation algorithm, SSA) can be approximated as a diffusion process, we 
may be interested in the dynamic features of the deterministic (drift) component of 
this diffusion. We perform traditional scientific computing tasks (integration, steady 
state and closed orbit computation, and stability analysis) on such a drift compo- 
nent using a SSA simulation of the Cyclic Lotka-Volterra system as our illustrative 
example. The results of short bursts of appropriately initialized SSA simulations are 
used to fit local diffusion models using A'ft-Sahalia's transition density expansions 
[1, 2, 3] in a maximum likelihood framework. These estimates are then coupled with 
standard numerical algorithms (such as Newton-Raphson or numerical integration 
routines) to help design subsequent SSA experiments. A brief discussion of the va- 
lidity of the local diffusion approximation of the SSA simulation (a jump process) 
is included. 



1 Introduction 

Reactive particle dynamic models arise in scientific fields ranging from phys- 
ical and chemical processes to systems biology [33, 34, 37, 41, 38, 19, 17]. 
Incorporating successive levels of detail in the modeling quickly leads to mod- 
els that are analytically intractable, necessitating computational exploration. 
Gillespie's Stochastic Simulation Algorithm (SSA) and its variants [41, 20, 19] 
have gained popularity in recent years for modeling so-called mixed reacting 
systems; the approach provides a middle ground between detailed molecular 
dynamics and lumped, Ordinary Differential Equation (ODE) descriptions of 
chemical kinetics, incorporating fluctuations. Knowing the kinetic scheme un- 
derlying such a simulation allows one to write, at the infinite particle limit, 
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the corresponding kinetic ODE. At intermediate particle numbers (N mol ), the 
SSA has been approximated with the continuous "chemical Langevin equa- 
tion" [19]. 

In what follows we will assume that the results of an SSA simulation can 
be successfully approximated through a continuous diffusion process. Explicit 
knowledge of the drift and noise components of such a process allows one to 
easily analyze certain features of the overall behavior; one might, for example, 
be interested in the bifurcation behavior of the "underlying" drift component 
of the model, including the number and stability of its steady states and their 
parametric dependence. In our work we assume that the only available simula- 
tion tool is a "black box" SSA simulator, in which the mechanistic rules have 
been correctly incorporated, but which we, as users, do not know: we can only 
observe the SSA simulator output. We want to perform a quantitative com- 
putational study of the underlying drift component. Since we cannot derive 
it in closed form (not knowing the evolution rules), we want to perform this 
study using the least possible simulation with the SSA code. The approach 
we use follows the so-called "equation- free" framework [28, 27]: in this frame- 
work traditional numerical algorithms become protocols for designing short 
bursts of numerical experiments with the SSA code. The quantities necessary 
for numerical computation with the unavailable model (time derivatives, the 
action of Jacobians) are estimated locally by processing the "fine scale" SSA 
simulations. In this work we extract such numerical information via paramet- 
ric local diffusion models using the transition density expansions proposed by 
A'it-Sahalia [2, 3]. The numerical procedures we illustrate can also be used, in 
principle, for different types of "fine scale" models if their output happens to 
be well approximated by diffusion processes. 

The article is organized as follows: In Section 2 we describe our illustrative 
model system. We then quickly outline the basic ideas underlying equation- 
free numerics (Section 3), and discuss our estimation procedure (Section 4). 
Our computational results are presented in Section 5, and we conclude with 
a discussion including goodness-of-fit issues. 

2 The Lattice Lotka-Volterra Model 

Our Cyclic Lotka-Volterra [36, 14] illustrative example consists of a three- 
species (X, Y and S) nonlinear kinetic scheme of the following form [36]: 

X + Y % 2Y 

Y + S H 2S 

S + X h 2X. 

In the remainder of the paper we will refer to it simply as LV. In the 
deterministic limit, this kinetic scheme gives rise to a set of three coupled 
nonlinear ODEs for the evolution of the concentrations X, Y and S. 
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d ^- = -k 1 XY + k 3 XS (1) 
dY 

— = -k 2 YS + kiYX 
at 

^ = ~k 3 SX + k 2 SY 
at 

The total concentration (X + Y + S) is constant over time; setting (without 
loss of generality) this constant to unity and eliminating 5* 

X + Y + S = l, => S = l-X-Y 

reduces the system to 



^ = X [k 3 - k 3 X - (h + k 3 )Y] (2) 
dY 

— = Y [-k 2 + (fcj + k 2 )X + k 2 Y] . 

For every (positive) value of k\,k 2 and £3 four fixed points exist: three 
trivial and one non-trivial steady state: 



(system invaded by S) 
(system invaded by X) 
(system invaded by Y) 
ki 

— (nontrivial fixed point) 
where 

K = fci + k 2 + k 3 . (3) 

An interesting feature of the phase space of the deterministic model is the 
existence of a one-parameter family of closed orbits surrounding a "center" 
(see Figure 3). The neutral stability of these orbits affects, as we will see 
below, the fixed point algorithms used to converge on them. The system is 
simulated through both the ODEs (2) and through an SSA implementation 
of the kinetic scheme (1) using k\,k 2 = 0.5 and k 3 = 0.7 throughout. 



X s = 0, Y s = 0, S = 1 
X s = 1, Y s = 0, S = 
X s = 0, Y s = 1, S = 

Y - -I V - tl G 



3 Equation Free Computation 

The basic premise underlying equation-free modeling and computation is that 
we have available a "black box" fine-scale dynamic simulator, and we believe 
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that an effective evolution equation exists (closes) for some set of (coarse- 
grained) outputs or observables of the fine scale simulation. As discussed in 
more detail in [27, 23], one can numerically solve this (unavailable explic- 
itly) equation through linking traditional numerical methods with the fine 
scale code; in particular, the classical continuum algorithms become proto- 
cols for the design of short, appropriately initialized numerical experiments 
with the fine scale code. The process starts by identifying the appropriate 
coarse-grained observables (sometimes also called order parameters); typically 
these variables are low-order moments of microscopically/stochastically evolv- 
ing particle distributions (e.g. concentrations for chemically reacting systems, 
like our example). In general, good coarse-grained observables are not known, 
and data analysis techniques to identify them from computational or experi- 
mental observations are the subject of intense current research [39, 8, 12]. If 
the unavailable "effective" equations are deterministic and reasonably smooth, 
short runs of the fine-scale simulator are used to estimate time derivatives of 
the coarse-grained observables; initializing fine scale simulations consistent 
with nearby values of the coarse-grained observables gives estimates of di- 
rectional derivatives (again assuming appropriate smoothness), and can be 
linked with matrix- free iterative linear algebra techniques (e.g. [26]). When 
an explicit evolution equation is available, these quantities, necessary in nu- 
merical computation, are obtained through function or Jacobian evaluations 
of the model formulas; here, they are estimated on demand from short com- 
putational experiments with the fine scale solver. If the underlying effective 
equation is stochastic, e.g. a diffusion, then the results of the short simulation 
bursts must be used to estimate both the drift and the noise components of 
the effective model - this is the case we study here. We will illustrate, using 
the SSA LV example, how certain types of computations can be accelerated 
by appealing to classical numerical methods. 

4 Estimation Procedure 

In what follows, we will assume that species concentrations are good observ- 
ables, and that the true process (the LV SSA simulation) can be adequately 
approximated by a diffusion process, that is, a stochastic differential equation 
(SDE) of the form: 

rfX t = n(X t :6)dt + £(X t ;0)dW t . (4) 

Here X t is a stochastic process which is meant to model the evolution of the 
observable(s), W t represents a vector of standard Brownian motions, and the 
functions /x(Xt;0) and £(Xt;#) are the drift and diffusion coefficients of the 
process. In the classical parametric setup, one assumes that the parameterized 
function families to which the drift and diffusion coefficient functions belong 
are known, and that the parameter vector 9 is finite dimensional. In practice 
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one rarely knows a class of functions which can be used to describe the global 
dynamics of the observables; in the equation free computations below, how- 
ever, we simulate the true process for only relatively short bursts of time. It 
therefore makes sense to (locally) consider the following SDE: 



where W t ,X t ,X , A and C £ M. d and B and D e R dxd (the d Brownian mo- 
tions are assumed independent; the vector multiplying them, by slight abuse of 
notation, contains the nonzero elements of the diagonal matrix S; extending 
to the correlated case is straightforward). 

This simple model is based on the fact that we expect smooth evolution 
of moments of the observables, while at the same time taking into account 
the state dependence of the noise (neglecting this dependence can cause bias 
in the estimation of the drift) . The parameters of this local linear model are 
estimated through techniques associated with maximum likelihood estimation 
(MLE). The motivation for using MLE techniques stems from the fact that 
under certain regularity conditions [42] such estimators are (asymptotically) 
efficient as regards the variance of the estimated parameter distribution. In 
addition, the asymptotic parameter distributions associated with MLE can 
sometimes be worked out analytically, or approximated through Monte Carlo 
simulations; this knowledge can guide the selection of the sample size necessary 
for a given desired accuracy in coarse-grained computations [11]. 

4.1 Maximum Likelihood Estimation for Discretely Observed 
Diffusions 

We now recall a few basic facts about MLE estimation; standard references 
include, e.g. [21, 24, 42]. It is assumed throughout that the exact distribu- 
tion associated with the parametric model admits a continuous density whose 
logarithm is well defined almost everywhere and is three times continuously 
diffcrcntiable with respect to the parameters [30] . 

MLE is based on maximizing the log-likelihood (£$) with respect to the 
parameter vector (for our model 6 = [A, B, C, D]): 



In the above equation, x corresponds to a matrix of observations £ ~R dxM 
where d is the dimension of the state and M is the length of the time series; 
/(x; 9) corresponds to the probability of making observation x. For a single 
sample path of a discretely observed diffusion known to be initialized at xo, 
/(x;#) can be evaluated as [21]: 





(6) 



M-\ 



f(x;6) = 6 X0 J] /(*ml 



x m _i;0). 



(7) 



m— 1 
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In this equation /(x TO |x TO _i;#) represents the conditional probability 
(transition density) of observing x m given the observation x TO _i for a given 
9 and 5 X0 is the Dirac distribution. In our applications, we search for the pa- 
rameter vector that is best over all observations (we have an ensemble of N 
paths of length M). In this case our expression for the log- likelihood (given 
the data and transition density) takes the form: 

N M 

^:=EE lo 4 /(x ™ |x ™- i;e) )- ( 8 ) 

i—l m— 1 

Assume the existence of an invertible symmetric positive definite "scaling 
matrix" matrix T^ M _ #) [31] associated with the estimator; the subscripts are 
used to make the dependence of the scaling matrix on M and 6 explicit. For 
the "standard" case N = 1 in time series analysis, under some additional 
regularity assumptions [24, 42], one has the following limit for a correctly 
specified parametric model: 

^ Mt§) (0 M -e)^N(O,I). (9) 

Here 8 is the true parameter vector; 8 M represents the parameters es- 

timatcd with a finite time scries of length M ; ==>■ denotes convergence in 
distribution [42, 21] under Pg (the distribution associated with the density 

/(x;#)); N(0, 1) denotes a normal distribution with mean zero and an iden- 
tity matrix for the covariance. For a correctly specified model family, T, M 
can be estimated in a variety of ways [44, 31]. The appeal of MLE lies in that, 
asymptotically in M, the variance of the estimated parameters is the smallest 
that can be achieved by an estimator that satisfies the assumed regularity 
conditions [24, 42]. 

4.2 Transition density expansions 

Here we briefly outline the key features of the recent work of A'it-Sahalia [2, 3] 
used in our coarse-grained computations below. The problem with using even 
a simple model like that given in equation 5 is that the transition density as- 
sociated with the process is not known in closed form. In recent years, many 
attempts to approximate the transition density have appeared in the liter- 
ature; some techniques depend on analytical approximations whereas others 
are simulation based (see, e.g. [1, 2, 5, 9, 18, 35]). We have used, with some 
success, the expansions found in [1, 3, 2]. High accuracy can be obtained us- 
ing this method to approximate the transition density associated with a scalar 
process; the multivariate case is discussed in [3]. The basic idea behind the 
scalar case, presented in [1, 2], is as follows: One first transforms the process 
given in equation 4 into a new process [2] : 
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J <r(u,t 



dY t = n Y (Y t ;6)dt + dW t (10) 
rX du 
1) 

^y(v;0)= (li) 
<7(7- 1 (2/;^);^) 2fe 17 (y ' U) ' U) 

An additional change of variables brings the transition density of the pro- 
cess closer to a standard normal density Z = A~?(Y — y ) where A is the 
time between observations. The transformations introduced allow the use of a 
Hermite basis set in order to approximate the transition density of the original 
process via the following series: 



p z (A,z\y o :0)^ (12) 

K 
J=0 



V v>(A,y ;6)= (13) 

OC 

i J Hj{z)pz(A,z\y ;e)dz:= 

— OO 

jp[Hi (AT* (Y t+A - yo )) \Y t = y ; 6} 

In the above, Hj represents the j th Hermite polynomial and is the 
standard normal density. The coefficients needed for the approximation are 
obtained through the conditional moments of the process Y t . A'it-Sahalia out- 
lines [2] a procedure which exploits the connection between the SDE and the 
associated Kolmogorov equations in order to develop a closed form expres- 
sion for the rfjp coefficients. The approximation is exact if K — > oo and the 
coefficient functions satisfy the assumptions laid out in [2]. In numerical ap- 
plications one must always deal with a finite K. Problems may arise in the 
truncated expansion: the approximation of the density may not normalize 
to unity or, worse, it may become negative (see [3, 1, 5] for some possible 
remedies). 

In the multivariate case, it becomes more difficult to introduce an analog 
of Yt [3]. Nonetheless, it is still possible to construct a series motivated by 
the methodology used in the scalar case; however, one now needs to expand 
in space and time, whereas the Hermite expansion yielded a series "in time 
only" [3]). A'it-Sahalia [3] outlines an approach which makes use of a recursion 
for calculating the coefficients of the expansion in the multivariate case. We 
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have had success in using these expansions, even in cases where convergence 
of the infinite series is not guaranteed by the conditions given in [1, 2]. Notice, 
for example, that our local models may allow a value of zero for the diffusion 
coefficient; using a different function class (made computationally feasible by 
the extension of Bakshi and Ju [5, 6]), such as sigmoidal functions for it, 
may help circumvent such problems. Other pathologies are discussed in [11]; 
estimates of the range of the parameters of interest [31, 42, 11] can enhance 
the algorithm performance. The comparison study [25] recommends the use of 
the expansions by A'it-Sahalia for a wide class of diffusion models. Beyond the 
estimation itself, these expansions can also be helpful in obtaining diagnostics 
that depend on knowledge of the transition density (such as goodness-of-fit 
tests [22]) and asymptotic error analysis [31]. 

5 Illustrations of Equation-free Computation 

Having estimated the parameters of a local model at a given state point opens 
the way to several computational possibilities. Such estimates, for example, 
can be used in an iterative search for zeroes of the (global, nonlinear) drift. 
A Newton-Raphson iteration for a (hopefully better) guess of this root in- 
volves the solution of set of linear equations for which both the matrix and 
the residual are available from the local linear drift. The resulting estimate of 
the root is then used to launch a new set of computational experiments with 
the "inner" SSA code, followed by a new estimation, linear equation solution, 
and so on to convergence. This illustrates the fundamental underpinnings 
of equation-free computation. Many numerical algorithms (here, root finding 
through Newton iteration) do not really require good closed-form global mod- 
els: each iteration only requires local information (the first very few terms of 
a Taylor series) in order to "design" the next iteration. Traditional contin- 
uum numerical methods can thus be thought of as protocols for the design of 
a sequence of model evaluations (possibly model and Jacobian evaluations, 
occasionally even Hessian evaluations). In the absence of an explicit formula 
for the model, the same protocol can be used to design appropriately ini- 
tialized computational experiments with a model of the system at a different 
level (here, the SSA simulator). Processing the results of these appropriately 
initialized short bursts estimate the quantities required for scientific computa- 
tion, as opposed to evaluating them from a closed-form model. The so-called 
"coarse projective integration" is another example of the same principle. Tra- 
ditional explicit integration routines require a call to a subroutine that eval- 
uates the time-derivative of a dynamic model at a particular state. In the 
absence of an explicit model, short bursts of simulations of a model of the 
system at a different level (again, here, SSA) can be used to estimate these 
time derivatives, and, through local linear models, extrapolate the state at a 
later time. The fundamental assumption underpinning this entire computa- 
tional framework, is that an explicit evolution equation exists, and closes, in 
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terms of the (known) coarse-grained observables of the fine-scale simulation 
(here, the concentrations of the SSA species). If this, unavailable in closed- 
form, equation is deterministic, then one only need to estimate a drift term 
from fine scale simulations; if, on the other hand, the coarse-grained equation 
is stochastic (fluctuations are important), then both the local drift and dif- 
fusion terms must be estimated. Certain computational tasks for stochastic 
effective, coarse-grained models require evaluations of both these terms (e.g. 
computations of stationary, equilibrium densities, or Kramers' type compu- 
tations of escape times for bistable systems, see for example [29, 23]). In this 
paper, we perform equation-free tasks for only the drift component of the 
model; sometimes it may be interesting to know whether the drift component 
dynamics possess zeros or closed loops, as well as their parametric depen- 
dence. Also, at infinite system size (practically, for sufficiently large particle 
numbers) the SSA actually closes as a deterministic ODE. 

Coarse Newton- Rap hson for the fixed point of the drift. 

In what follows we work at system sizes large enough that a diffusion approxi- 
mation of the SSA output is meaningful, and -even more- the dynamics of the 
drift component of the diffusion are close to the kinetic ODE scheme dynam- 
ics. The neutral stability of the fixed point and the closed loops of the kinetic 
ODE suggest comparable features for the estimated drift, which we set out 
to investigate. We find the nontrivial root of the estimated drift F(X;0) = 1 
through a coarse Newton- Raphson procedure as follows: An ensemble of N pat h 
SSA simulations are initialized in a neighborhood of the current guess X of 
the root. Each is evolved in time, and the simulations are sampled uniformly 
M times during a time interval of length r. A local SDE model of the type 
(5) is estimated using the transition density expansions of Ai't-Sahalia in an 
MLE-type scheme; the resulting model parameters are used to update the 
root guess through 

X„ = X„_i ' |x=x n _iF(X„_i; 6) « X„_i — B 1 A. (14) 

Figure 1 shows this procedure for two different values of N pat h (other 
parameters are noted in the caption). Newton-Raphson type procedures for 
isolated roots are known to converge quickly given a good initial guess; fur- 
thermore, upon convergence, the eigenvalues of the linearization of the drift 
are contained in the matrix B. Estimates of these eigenvalues for different 
Npath are listed, upon convergence of the root finding procedure, in Table 1. 
The equation- free iterates approach the deterministic ODE root (see inset); 
the latter is known to possess two pure imaginary eigenvalues. The estimated 
(from local models) eigenvalues are also characterized by a relatively small 
(0(1CT 2 )) real part. 

1 We use F(-; 9) to denote the right hand side of a general deterministic ODE; here 
F(-;#) is the estimated n(-;8) 
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-11 1 1 1 ' 

1 2 3 4 5 

Iteration 

Fig. 1. Coarse NR to find stationary points. Roots of F(X) = %r- using estimated 
(local) linear SDEs. Parameters: iV mo! = 1 x 10 4 , r = 5.032928126 x 10" 1 , M = 300. 
N pa th values are shown in the legend and the I 2 distance of the current guess from 
the deterministic ODE root is shown in the inset. Z\X* (inset y-axis) represents the 
difference between the current guess and the steady state of the ODE. 

Table 1. Representative real and imaginary parts of the eigenvalues of the estimated 
drift upon convergence to the nontrivial fixed point X w (0.2941,0.41176); the 
deterministic ODE solution has a pair of pure imaginary eignenvalues. 





Re 


Im 


N path = 100 


-2.53 x 10~ 2 


3.01 x 10 _1 


N path = 600 


-2.39 x 1Q-' J 


2.63 x 10" 1 



Coarse Projective Integration for the drift 

A variety of numerical integration algorithms can be implemented in our 
framework. Single step methods of the general form 



X n = X n _ 1 + *(X n _ 1 ,X„;^t). 



(15) 
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include the explicit and implicit Euler algorithms, (for which 3? is Z\tF(X n _ 1 ) 
and AtF(K ) respectively). Estimates of the drift at X can be immediately 
used in a "coarse forward Euler", while the estimated B can be used in a 
root-finding procedure, along the lines illustrated above, in a "coarse back- 
ward Euler" scheme. Other schemes can be simply implemented. Here we 
only demonstrate (explicit) coarse forward Euler; predictor-corrector schemes 
(more appropriate for stiff problems) are illustrated in [10]. Representative re- 
sults for our LV problem are shown in Figure 2. The deterministic ODE trajec- 
tory (dashed lines connecting points) is compared to the projective integration 
of the drift component of an SDE estimated locally from SSA simulation en- 
sembles. One clearly sees the evolution of the ensemble of SSA trajectories 
initialized at every numerical integration point; N pat h such trajectories were 
evolved and observed uniformly M times over a time interval r. The results 
were processed through the estimation scheme and the value of the drift at 
the original point Xo provided the forward Euler estimate of the "next" point 
through Xi = X + Z\tF(X ). The procedure is then repeated. 

Several algorithmic parameters must be carefully selected in such compu- 
tations. In our case the "lifting" problem (the initialization of SSA simulations 
at a given value of the coarse observables) is straightforward because of the 
"mixed" nature of the SSA simulation; in general, the successful initialization 
of a fine scale code consistent with a few coarse observables can be a com- 
plicated and difficult issue, requiring, for example, preparatory constrained 
dynamic runs [4, 40]. 

Another important parameter is the length of the integrator "projective" 
step, At, which for deterministic problems is set by stability and accuracy 
considerations. Stability discussions for projective integration can be found in 
[28] ; here the issue is complicated by the fact that the model is estimated rather 
than evaluated. Multiscale methods for SDEs, including error estimates, can 
be found in the work of [43, 16]. The total "microscopic integration time" 
denoted by r and the time between observations = 5t := also require 
careful selection. If r is too large, the simple linear model may break down as 
nonlincaritics in the real system manifest themselves. If the assumed diffusion 
model is correct, there is no upper limit on M; yet a diffusion approximation 
of a different underlying process, such as the jump SSA here, will break down 
if the data is sampled too frequently. Similar issues have been addressed in the 
control literature [15]. Later on we will outline a goodness-of-fit test that can 
be used to guide the selection of such algorithmic parameters. In this work, 
short SSA trajectories in each ensemble are initialized at the same base point 
Xo, or uniformly in a small neighborhood around it; we have not yet explored 
optimal initialization. 

Equation-Free Coarse Variational Calculations 

A slight extension of the above coarse integration procedure is the implemen- 
tation of equation free integration of variational equations. The need for these 
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0.45 



Fig. 2. Illustration of Coarse Projective Integration. N pat h = 600, N mo i = 1 x 10 4 , 
At = .50329, r = ^. 



arises naturally in our example when we attempt to construct an algorithm 
that searches for possible closed orbits in the dynamics of the estimated drift, 
and attempts to converge on them. Closed orbits that are limit cycles can 
be found as (isolated) fixed points of an appropriate Poincare map. In the 
deterministic LV problem, however, one has a one-parameter family of such 
orbits, and the fixed points of the Poincare map are not isolated. Anticipating 
a family of such closed orbits for our estimated drift model, we isolate a sin- 
gle orbit from this one-parameter family by selecting its period (the Poincare 
return time). 

For a deterministic model, the initial value problem for the variational 
equations is 



F(X;0) 



dX 
~~dt 

X(t = 0) = X IC 
dV _ <9F(X;0) 

~dt ~ ax ' 



(16) 



V 



Coarse-graining SSA LLV models via local diffusion models 13 
V(t = 0) = I. 

If X £ R d then V £ R dxd . We use the results of integrating such variational 
equations to locate closed orbits as zeroes of the equation G(X) = X 3> T (X) 
where 3?r(') represents the result of integration from the (deterministic) initial 
condition X /<7 for time r. To isolate the zeroes we seek, we select a Poincare 
plane through the value Xp = 0.3 of the first coordinate, and the return 
time; we thus have one equation with one unknown, the Y coordinate of the 
intersection of our particular closed orbit with the chosen Poincare plane. For 
our coarse integration, the return time r is typically too large to permit a 
single local diffusion model to accurately describe the dynamics; we therefore 
use the following procedure: 

• Specify r and the number N gri d of local models we will use along the orbit, 
each valid for TY lacro := 



• Simulate N pat h SSA trajectories starting at the current fixed point guess; 
use the data as above to estimate the first local linear model. Use its drift 
(and the matrix B) to obtain the next "base point" as well as to step the 
variational equations for time 7y iacro . 

• Repeat N gri d times (see Figure 3). 

The output of this procedure gives us the residual of the fixed point equa- 
tion we wish to solve; the results of the variational integration at time r 
(which, upon convergence, will give us an estimate of the monodromy ma- 
trix) are then used to compute the Jacobian of the fixed point scheme. One 
Newton-Raphson step for the Y coordinate of the fixed point is taken, and 
the procedure is then repeated. Representative numerical results are shown in 
Figure 4. Because of the neutral dynamics, the eigenvalues of the monodromy 
matrix upon convergence are both equal to 1 (in the deterministic ODE). Ta- 
ble 2 shows representative eigenvalue upon convergence for different N pa th 
(sometimes the eigenvalues are numerically found as complex conjugates with 
a small imaginary part). Clearly, in addition to the algorithmic parameters 
involved in coarse projective integration, we should now also take into account 
the desired accuracy of the variational integration (quantified in part by the 
existence of an eigenvalue equal to unity upon convergence). 



Table 2. Representative monodromy matrix eigenvalues upon convergence of the 
fixed point iteration for two distinct N pat h computations (see text) . 



N path = 100 


(0.9698,0.1707) 


(0.9698, -.1707) 


N path = 400 


(0.9159,0.0252) 


(0.9159, -.0252) 
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Fig. 3. Coarse closed orbit computations for the Lotka Volterra model. The de- 
terministic model phase portrait contains an infinite number of closed orbits. Three 
such deterministic orbits (obtained by Runge-Kutta integration) are plotted here. To 
find the closed orbit with a specified period r, we use the Poincare surface Xp — 0.3, 
shown as a solid line. The Jacobian of the coarse Newton-Raphson scheme is com- 
puted through variational integrations based on the estimated drift from ensembles 
of SSA simulations initialized at the N gr id base points shown (see text). 

6 Discussion 

We have illustrated the implementation of certain coarse-grained computa- 
tions with the LV model; one focus was the coarse-grained integration of 
the variational equations for the SSA-based drift estimation, as well as the 
modifications of the coarse Newton-Raphson iteration dictated by the neutral 
stability of the dynamics (the existence of infinitely many closed orbits in the 
ODE limit, which appears to approximately persist in our computations). The 
second focus was the use of A'it-Sahalia's expansions to estimate local linear 
SDEs from short bursts of SSA data as an intermediate step. This naturally 
leads to some crucial questions about the goodness-of-fit of the simple SDE 
models: (a) is the diffusion approximation a "good" description of the dynam- 
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Fig. 4. Coarse Newton-Raphson for finding closed orbits of a specified period. The 
zeroes of G(Y) = Y — & T (Y) were calculated using a Jacobian evaluated from coarse 
variational integration based on SSA simulations. Parameters: N mo i — 1 x 10 4 , 
t = 2.0131712504 x 10\M = 300, N grid = 40 (N pat h given in the legend). The 
initial guess was Y = 0.53. For the deterministic ODE model the fixed point is 
Y ODE w 0.518; the coarse fixed point for N path = 400 was calculated to be Y SSA « 
0.5075. 



ics? (b) Is the linear approximation valid for the time series length chosen? 
and (c) How reliable is the model for making predictions/forecasts ? 

One should quantitatively know how large N mo i needs to be, for a given 
sampling frequency, for a diffusion model to be a statistically meaningful ap- 
proximation [7, 19]. Sampling too often may be detrimental in many diffusion 
approximations (e.g. [15]). Local linear models (i.e. short truncations of Tay- 
lor series) are used extensively in scientific computations, but only for short 
time steps, whose length is determined by overall stability and accuracy con- 
siderations. Similar considerations arise in choosing the r used for SSA data 
collection towards the estimation of the local linear SDE models used here; 
clearly, when the underlying drift is nonlinear, r cannot be too large. A useful 
diagnostic tool for questions (a) and (b) applicable if one does have an ac- 
curate transition density approximation, is the probability integral transform 
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Fig. 5. Evolution of actual and model densities. The top figure shows the evolution 
of the SSA process, initialized as a Dirac distribution (marked in red); the bot- 
tom plot shows an ensemble of numerical simulations of the ideal diffusion model 
using the parameters estimated from the SSA. Both distributions are plotted at 
j, §, ^,and r. Relevant parameters: N pat h = 5 x 10 3 , M = 300, N mo i = 1 x 10 4 , 
r = 0.50329. 

[13, 22]. Using the data and the (assumed known) exact transition density, 
one creates a new random variable which, for a correctly specified model, 
has a known distribution. The method is applicable to both stationary and 
non-stationary time series; furthermore it depends on integrations of the tran- 
sition density approximation rather than differentiations. Given the data, one 
(appropriately) estimates model parameters and then constructs the random 
variables Z n for each observation 2 (x n ). The construction below follows that 
in Section 3 of [13]: 

Z n ■= J p(x' n \x n -i;0)dx' n 

— oo 

2 The method applies to both a vector and scalar process, however the construction 
is easiest to demonstrate in the latter case. See [22] for the multivariate extension. 
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Z n ~ q(Z n ) 
Xn ^ f(%n\3'n— l) 



d¥(x n \x n - 1 ) 



dx r 



The symbol ~ denotes that the random variable on the left of the symbol is 
distributed according to the density to the right. Under a correctly specified 
model, the Z n 's are independent and uniformly distributed on [0, 1], indepen- 
dent of the transition density [13]. In [22] a comprehensive suite of statistical 
tests are reviewed which exploit knowledge of the transition density and the 
transformation shown above. Figure 6 plots a kernel density estimate (see 
equation 6 on page 44 in [22]) which is based on the estimated parameters 
and the observed data. If the model is correctly specified, the infinite sample 
size density should be the product of two uniform densities. Test statistics can 
be created from this function (see [22] for details). 

Inspection of the figures shows that, for a particular representative SSA 
ensemble run for N mo i = 1 x 10 4 , and a particular sampling frequency, the 
diffusion approximation is not acceptable; the situation appears better for 
Nmoi > 4 x 10 5 . It is interesting to notice that, while N mo [ = 1 x 10 4 is not 
large enough for the conditions of Figure 5, visual inspection of the empirical 
and the SDE-based density evolution might suggest otherwise. In traditional, 
continuum numerical algorithms issues of on-line error estimation, time-step 
and mesh adaptation are often built-in in modern, validated software. There 
is a clear necessity for incorporating, in the same spirit, hypothesis testing 
techniques in codes implementing the type of computations we described here; 
yet automating such processes appears to be a major challenge. 

In our next application, we evolved an ensemble of trajectories starting 
from a Dirac initial distribution, and then recorded the Poincare map for 
each individual trajectory over a long simulation period. Figure 7 shows the 
evolution of the Y coordinate of these trajectories as function of the map 
iterate. For long times, different initial conditions in the ensemble approach 
some of the "extinction" fixed points of the ODE vector field (see the vertical 
lines in Fig. 7); once there, the system no longer changes over time. Visual 
inspection of the evolution of the ensemble suggests that one might try to 
coarse-grain the Poincare map evolution as a model SDE; the insets in the 
figure show the initial evolution of the mean and the variance of the Poincare 
map iterates. The smooth line in the insets, a simple least squares fit, seems 
to suggest a systematic evolution towards "larger" oscillations, bringing the 
system closer to extinction. If this evolution could be well approximated locally 
by a diffusion processes, approximations similar in spirit to the ones shown in 
this article might be used to explore features of the distribution of extinction 
times for the problem. 

Acknowledgements This work was partially supported by a Ford Foun- 
dation/NRC Fellowship to (CC) and an NSF ITR grant (IK). 



18 C. P. Calderon, G. A. Tsekouras, A. Provata, and I. G. Kevrekidis 




Ul 04 03 Cfl. O.T Oft O.B 



Fig. 6. Towards hypothesis testing. The function plotted corresponds to an empiri- 
cal estimate of the two-dimensional density function described in [13, 22]. The data 
are obtained from the same ensemble of SSA simulations as in Fig 5; the top figure 
is for N mo i — 1 x 10 4 and the bottom for N mo i = 4 x 10 5 . In the infinite sample 
limit and for a correctly specified model the density would be unity in the entire 
support ([0, 1] x [0, 1]) of the function. The figure suggests that the observations of 
the larger system are closer to a diffusion model. 
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